Problème torsion d’arbres

Librairies

Show code cell content Hide code cell content

# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue

Géométrie

Définition de la géométrie et du maillage

Show code cell content Hide code cell content

# Initialisation de GMSH
gmsh.initialize()
gmsh.model.add("ellipse")

# Définition des paramètres géométriques et de maillage
a = 0.2  # Demi-grand axe de l'ellipse en mètres
b = 0.2  # Demi-petit axe de l'ellipse en mètres
h = 1.0  # Hauteur du cylindre elliptique en mètres
# Calcul du rayon équivalent du cylindre elliptique
R = 1/2 * np.sqrt(a**2 + b**2)

# Création de la géométrie : ellipse de base
ellipse = gmsh.model.occ.addEllipse(0, 0, 0, a, b)

# Création d'une boucle fermée à partir de l'ellipse
curve_loop = gmsh.model.occ.addCurveLoop([ellipse])

# Création d'une surface plane à partir de la boucle
surface = gmsh.model.occ.addPlaneSurface([curve_loop])

# Extrusion de la surface pour créer un cylindre elliptique
gmsh.model.occ.extrude([(2, surface)], 0, 0, h)

# Synchronisation du modèle géométrique
gmsh.model.occ.synchronize()

# Configuration des paramètres de maillage
lc = 0.02  # Taille caractéristique des éléments du maillage en mètres
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)

# Génération du maillage volumique (3D)
gmsh.model.mesh.generate(3)

# Sauvegarde du maillage au format .msh
gmsh.write("ellipse.msh")

# Finalisation de GMSH
gmsh.finalize()

# Lecture du fichier .msh avec meshio
msh = meshio.read("ellipse.msh")

# Conversion et sauvegarde du maillage au format XDMF (compatible avec FEniCSx)
meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))

# Lecture du maillage XDMF avec FEniCSx
with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
    domain = xdmf_file.read_mesh(name="Grid")
    domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

# Affichage des informations sur le maillage
print("Maillage créé et importé avec succès.")
print(f"Nombre de sommets : {domain.topology.index_map(domain.topology.dim).size_global}")
print(f"Nombre d'éléments : {domain.topology.index_map(domain.topology.dim-1).size_global}")

Visualisation de la géométrie

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=700, height=400)
interactive_panel
from myst_nb import glue
glue("img", interactive_panel,display=False)

Configuration du problème

import ipywidgets as widgets
from IPython.display import display

# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)

# Code Fenics ou autre qui dépend de "rayon.value"

import ipywidgets as widgets
from IPython.display import display

# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)

# Code Fenics ou autre qui dépend de "rayon.value"

Définition de l’espace fonctionnel

# Exemple de dropdown pour choisir un matériau
materiau = widgets.Dropdown(
    options=['Acier', 'Aluminium', 'Plastique'],
    value='Acier',
    description='Matériau:',
)

# Fonction pour mettre à jour le slider en fonction du matériau choisi
def update_slider(change):
    if change['new'] == 'Acier':
        rayon.value = 0.05  # Valeur par défaut pour l'acier
    elif change['new'] == 'Aluminium':
        rayon.value = 0.04  # Valeur par défaut pour l'aluminium
    else:
        rayon.value = 0.03  # Valeur par défaut pour le plastique

materiau.observe(update_slider, names='value')
imk=materiau.observe(update_slider, names='value')

display(materiau)

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
from myst_nb import glue
glue("img", interactive_panel)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[101], line 3
      1 # Utiliser Panel pour créer un rendu interactif
      2 pn.extension('vtk')
----> 3 interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
      4 from myst_nb import glue
      5 glue("img", interactive_panel)

AttributeError: 'NoneType' object has no attribute 'ren_win'

Définition des frontière du domaine

Show code cell content Hide code cell content

# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
    """Identifie la face inférieure du cylindre (z = 0)"""
    return np.isclose(x[2], 0)

def top_face(x):
    """Identifie la face supérieure du cylindre (z = h)"""
    return np.isclose(x[2], h)

def lateral_face(x):
    """
    Identifie la surface latérale du cylindre elliptique
    Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
    """
    tolerance = 1e-5  # Tolérance pour la comparaison numérique
    return (np.isclose((x[0]**2 / a**2 + x[1]**2 / b**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)

# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1

# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)

# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_h])  # Utilisation uniquement de la face supérieure
marked_values = np.hstack([np.full_like(Sigma_h, 1)])  # Marqueur 1 pour la face supérieure

# Tri et création des tags de maillage
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, domain.topology.dim - 1, 
                                  marked_facets[sorted_facets], 
                                  marked_values[sorted_facets])

# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)

Visualisation des frontière du domaine

Show code cell outputs Hide code cell outputs

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc1 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc2 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc3 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau

# Définition des constantes élastiques du matériau
mu = 80e9      # Module de cisaillement (G) en Pa
lambda_ = 120e9  # Premier paramètre de Lamé (λ) en Pa

Définition des chargements

α = 20.0 m^-1
C = 1005309649.1487346 kg.m^2/s

Résolution

Show code cell content Hide code cell content

# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
    """Tenseur de déformation linéarisé"""
    return ufl.sym(ufl.grad(u))

def sigma(u):
    """Tenseur des contraintes (loi de Hooke)"""
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx  # Forme bilinéaire
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(1)  # Forme linéaire

# Définition et résolution du problème linéaire
problem = LinearProblem(a, L, bcs=[bc1], 
                        petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
uh = problem.solve()

Visualisation des résultats

Visualisation des déplacements

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales

Visualisation des déplacements amplifiés

Visualisation des déplacements normalisés et selon les axes principaux

Visualisation des composantes du tenseur des déformations

Visualisation des composantes du tenseur des contraintes

Calcul des contraintes de von Mises

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats

Calcul des valeurs maximales

Show code cell content Hide code cell content

# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))

# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)

# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)

# Calcul de l'angle de torsion théorique
theoretical_twist = ang * h  # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre

# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100

# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100

Moment d'inertie polaire J : 6.283185e-04 m^4
Angle de torsion théorique θ : 20.000000 radians
Angle de torsion théorique θ : 1145.92 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 19.888688 radians
Rotation Z maximale simulée : 1139.54 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 19.901232 radians
Norme de rotation maximale : 1140.26 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 124.69 degrés
RX max  : 2.176298 radians
RY max : 122.15 degrés
RY max  : 2.131950 radians